Comparison of computer simulations and clinical treatment results of magnetic resonance‐guided focused ultrasound surgery (MRgFUS) of uterine fibroids

Abstract Purpose Magnetic resonance‐guided focused ultrasound surgery (MRgFUS) can be used to noninvasively treat symptomatic uterine fibroids by heating with focused ultrasound sonications while monitoring the temperature with magnetic resonance (MR) thermometry. While prior studies have compared focused ultrasound simulations to clinical results, studies involving uterine fibroids remain scarce. In our study, we perform such a comparison to assess the suitability of simulations for treatment planning. Methods Sonications (N = 67) were simulated retrospectively using acoustic and thermal models based on the Rayleigh integral and Pennes bioheat equation followed by MR‐thermometry simulation in seven patients who underwent MRgFUS treatment for uterine fibroids. The spatial accuracy of simulated focus location was assessed by evaluating displacements of the centers of mass of the thermal dose distributions between simulated and treatment MR thermometry slices. Temperature–time curves and sizes of 240 equivalent minutes at 43°C (240EM43) volumes between treatment and simulation were compared. Results The simulated focus location showed errors of 2.7 ± 4.1, −0.7 ± 2.0, and 1.3 ± 1.2 mm (mean ± SD) in the anterior–posterior, foot–head, and right–left directions for a fibroid absorption coefficient of 4.9 Np m–1 MHz–1 and perfusion parameter of 1.89 kg m–3 s–1. Linear regression of 240EM43 volumes of 67 sonications of patient treatments and simulations utilizing these parameters yielded a slope of 1.04 and a correlation coefficient of 0.54. The temperature rise ratio of simulation to treatment near the end of sonication was 0.47 ± 0.22, 1.28 ± 0.60, and 1.49 ± 0.71 for 66 sonications simulated utilizing fibroid absorption coefficient of 1.2, 4.9, and 8.6 Np m–1 MHz–1, respectively, and the aforementioned perfusion value. The impact of perfusion on peak temperature rise is minimal between 1.89 and 10 kg m–3 s–1, but became more substantial when utilizing a value of 100 kg m–3 s–1. Conclusions The results of this study suggest that perfusion, while in some cases having a substantial impact on thermal dose volumes, has less impact than ultrasound absorption for predicting peak temperature elevation at least when using perfusion parameter values up to 10 kg m–3 s–1 for this particular array geometry, frequencies, and tissue target which is good for clinicians to be aware of. The results suggest that simulations show promise in treatment planning, particularly in terms of spatial accuracy. However, in order to use simulations to predict temperature rise due to a sonication, knowledge of the patient‐specific tissue parameters, in particular the absorption coefficient is important. Currently, spatially varying patient‐specific tissue parameter values are not available during treatment, so simulations can only be used for planning purposes to estimate sonication performance on average.


Implementation details of a simulated sonication
Calculation of absorbed power density The surface velocities on the transducer elements were initially set to have magnitude 1, and the phases of the velocities on the transducer were calculated using a speed of sound value of 1540 m/s found in the treatment planning system configuration file, by back propagating from a target point calculated by using treatment parameters of the treatment planning system log and configuration files.
To take into account the acoustic power used to drive the transducer, the resulting absorbed power density fields were scaled based on the nominal power value in the treatment planning system log files in such a way that the power emitted by the transducer surface calculated via matched the nominal power value in treatment planning system log files which is assumed to be approximately equal to brush target measurements i.e. the absorbed power density fields calculated using a velocity of 1 were multiplied by . It was assumed that the power calibration of the transducer was performed while the transducer was immersed in oil, and thus the attenuation of the oil was modeled as zero in the acoustic simulations. The sizes of the areas of the contours on which the velocity was calculated were chosen to be large enough that the surface normal vectors of the transducer elements intersected the contours. The spatial discretization was picked to be 5 points per mm for each of the contours. This discretization value corresponds to approximately 6-7 points per wavelength for the 1.2 MHz sonications and 5-6 points per wavelength for the 1.44 MHz sonications. These more finely discretized contours were obtained by interpolating and then smoothing the contours manually segmented from the MR images. The skin-gelpad and gelpad-oil contours were derived from the smoothed skin-fat interface as described in the section titled 'Simulation Geometry, Patient Anatomy, and Coregistration'. In the phantom simulations, the phantom-membrane contour separating the phantom from the oil that the transducer is immersed in was manually defined/determined, followed by coregistration, interpolation and smoothing as was described for the manually segmented anatomical contours for the patient simulations.
The pressure simulation calculation volume was chosen to have dimension of approximately 7 by 7 by 9 cm and was centered about the treatment cell position with the misregistration correction value added to it. To save computational time, the pressure field was calculated using a multi discretization grid technique (Ellens and Hynynen 2014), beginning with a spatial discretization length of 0.8 mm, which was narrowed to 0.4 mm, and finally to 0.2 mm by interpolating and recalculating interpolated values in regions having a high Laplacian value.
For each trajectory, the absorbed power density field of a target point is first calculated using the targeting parameters of that point on the trajectory. The absorbed power density fields of the remaining points on the trajectory were obtained by rotation. The point of rotation was taken to be the misregistration corrected cell position with an additional lateral shift equal to the amount of lateral displacement of the calculated absorbed power density maximum from its corresponding target point. The rotation of a field calculated within a rectangular volume geometrically introduces points which are not well defined around the boundaries within the periphery of that volume, which were set to zero in all of the absorbed power density fields for the trajectory, resulting in fields trimmed to a cylinder like shape.

Thermal simulations
The temporal discretization dt was set to 12.5 ms and the spatial discretization used was 0.2 mm. The temperature simulation volume was chosen to be coterminous with the ultrasound volume, and was padded with additional voxels that were modeled to have constant temperature. The padding was removed after the simulation. An initial temperature condition of 37 degrees Celsius was assumed for the entire volume. A constant temperature boundary condition was implemented by holding the outermost two voxels on each face of the padded volume at fixed temperature. On each timestep, for each voxel that was not held at constant temperature, the finite difference calculation used two spatial neighbors on both sides of the voxel in each of the three grid directions.
Whichever media from the layered model were within the thermal simulation volume, were taken into account using tissue parameters of table 1. The thermal dose was calculated at each time step and perfusion was set to zero once the dose reached 240 equivalent minutes at 43 °C. Temperature rise was not simulated in oil meaning that any voxels consisting of oil were kept at constant temperature.
Trajectory point switching time of 50 ms, and the order of sonication of points on the trajectory used were derived from the treatment planning system configuration files. The durations of sonication trajectories were approximated by parameters in the treatment planning system log files and clinical treatment reports. A cooling simulation was performed having duration approximately from the end of sonication to the end of scanning that was derived from timestamps within the log files.

MR-thermometry simulation
The MR-thermometry simulation was done by first converting the simulated temperature values in the simulation grid having 0.2 mm discretization to phases via where is the temperature sensitivity coefficient, is the gyromagnetic ratio, is the magnetic field strength and TE is the echo time. Next, a temperature value for each MR-thermometry voxel was obtained by calculating the temperature using the same equation, but from the phase value of the average of the complex proton signal within that voxel. The numerical values of ppm °C -1 , MHz Tesla -1 , Tesla, and ms were used. For each sonication, a coronal stack and sagittal slice were placed at the treatment cell position. Another analysis was also carried out in which a coronal stack and sagittal slice were placed at the treatment cell position with the misregistration corrected value added to it. Unless otherwise stated, results shown are performed with the coronal stack and sagittal slice placed at the treatment cell position. The word 'placed' in this context means the following: In the case of the coronal stack 'placed at the treatment cell position' means that the center of the central voxel of the second coronal slice is placed at the treatment cell position. In the case of the sagittal slice 'placed at the treatment cell position' means that the voxel grid of the sagittal slice is defined in such a way that a voxel of the slice has its center at the treatment cell position (this voxel might not be the central voxel of the slice in the sense that the spatial extent of the slice in the anterior posterior direction with respect to the treatment cell position may be asymmetric due to the size of the simulation volume). The sampling interval of the simulated MR-slices was set to 2.57 s to coincide with the sampling interval of slices acquired during treatments. The voxel dimensions of the slices were set to 2.083 by 2.083 by 7 mm to match the ones used in the clinical slices.

Comparison of Clinical Data to Simulation Results:
To assess the accuracy of the simulations pertaining to magnitude of temperature rise, two metrics were used: peak temperature curves as a function of time evaluated using the sagittal slice of simulation and of treatment as well as thermal dose threshold volume sizes. The 240EM43 focal volume thermal dose distribution sizes were compared by comparing the 240EM43 thermal dose focal volumes obtained via voxel summation of clinical and simulated coronal slices. The overlapping information in the sagittal slice is not used in the evaluation of thermal dose: There is only one sagittal slice (7mm thick) which is not thick enough to cover a treatment cell having a 12mm diameter, furthermore the voxels (having already undergone spatial averaging) of the sagittal slice are oriented perpendicular (7mm in LR direction) to those of the coronal slices (7mm in AP direction). When presenting data involving averages of many sonications, the data is presented in relative quantities for purposes including power normalization as well as ease of presentation. However, several examples involving these metrics are presented using absolute quantities as they have the advantage of providing additional information thus enabling presentation of more informative concrete examples of clinical relevance (e.g. peak temperature curve).
Temporal registration of clinical and simulated MR-thermometry slices Peak temperature time curves calculated from the clinical and simulated sagittal MRthermometry slices were compared. The first simulated MR-thermometry slice was simulated to occur at approximately 2.57 seconds after the beginning of sonication, whereas the first sagittal acquisition of a sonication of a patient treatment occurred at some time before 2.57 seconds after sonication. The sampling interval between slices was approximately 2.57 seconds for both simulated and clinical slices. However, as a practical matter, on the figures where temperature is represented as a function of time, the zero value of the time axis is set to the beginning of the simulated sonication and the treatment sagittal dynamics are temporally registered (shifted in time) in such a way that the first temperature acquisition is set to occur at t=0, i.e. The first treatment temperature acquisition dynamic is set to have a time value of 0 seconds so that simulated and treatment temperature curves of any given sonication have the same set of time coordinates.

Verification of Transducer Element Phase Calculation
To verify that the transducer element initial velocity phases were appropriately calculated via backpropagation using the targeting parameters, pressure field simulations were ran in homogenous media with initial velocities of transducer elements having the negatives of the transducer element velocity phases in generator reports printed using the Sonalleve system software during the phantom experiment and with initial velocities calculated using the targeting parameters for three sonications. Sonications with three different values of electronic steering were simulated. These sonications were performed at three different depths. The treatment cell positions were selected to be approximately 1 cm apart from one another in depth. For these sonications, the transducer position was approximately the same (varied by less than 1 mm in each of three directions), and the electronic steering value in the depth direction used in the targeting parameters ranged from 4 mm in the anterior direction to 20.9 mm in the posterior direction. When the results of each simulated sonication using the two aforementioned initial velocities were compared, in each case the locations of the absorbed power density maxima coincided, the magnitudes of the absorbed power density maxima differed by less than 1 % and the linear absorbed power density profiles through the absorbed power density maximum along the three axes were in good agreement (see supplementary figure S1) suggesting that the backpropagation calculation using the targeting parameters generated transducer element phases that yield practically the same acoustic fields as the phases used to drive the transducer during the experiment.

Effect of variation of perfusion on center of mass of thermal dose distribution
When comparing the anterior-posterior location of the center of mass of a thermal dose distribution of the sagittal slice of a sonication simulated using a value of perfusion of 1.89 kg(m -3 s -1 ) to that of a higher value the results are the following: regarding the N = 20 sonications simulated using a fibroid absorption value of 8.6 Np(m -1 MHz -1 ) in which perfusion was varied to take on the values of 1.89, 10 and 100 kg(m -3 s -1 ), changing omega from 1.89 to 10 kg(m -3 s -1 ) yielded displacements ranging from 0 to 0.3 mm towards the posterior direction, and changing omega from 1.89 to 100 kg(m -3 s -1 ) yielded displacements ranging from 0.7 mm towards the anterior direction to 1.4 mm towards the posterior direction; regarding the N=2 sonications simulated using a fibroid absorption value of 1.2 Np(m -1 MHz -1 ) in which perfusion was varied to take on the values of 1.89 and 10 kg(m -3 s -1 ), changing omega from 1.89 to 10 kg(m -3 s -1 ) yielded displacements of less than 0.1 mm; regarding the N=5 sonications simulated using a fibroid absorption value of 1.2 Np(m -1 MHz -1 ) in which perfusion was varied to take on the values of 1.89, 10 and 100 kg(m -3 s -1 ), changing omega from 1.89 to 100 kg(m -3 s -1 ) yielded displacements ranging from 8.3 mm to 44.8 mm towards the anterior direction and changing omega from 1.89 to 10 kg(m -3 s -1 ) yielded displacements ranging from 0.2 mm to 6.1 mm towards the anterior direction.
In some sonications, including some of the aforementioned sonications for which perfusion was varied, simulations utilizing the fibroid absorption value of 1.2 Np(m -1 MHz -1 ) had a thermal dose center of mass on the sagittal slice that based on numerical value appears to be displaced further anteriorly from the non-simulated thermal dose center of mass than those of the simulations utilising larger values of absorption. This apparent displacement at least in the case of some of the aforementioned sonications for which perfusion was varied, is likely not due only to properties of the thermal dose distribution of the focal volume, but likely instead is caused by a lack of heating to such an extent in the focal volume that the thermal dose distribution of the nearfield of the simulated sagittal slice largely contributes to the apparent displacement. The same phenomenon due to lack of focal heating may apply for simulations utilizing a fibroid absorption value of 18 Np(m -1 MHz -1 ) as well. Visual inspection of the last simulated sagittal thermal dose distribution slice of the five sonications simulated using a fibroid absorption value of 1.2 Np(m -1 MHz -1 ) in which perfusion was varied to take on the values of 1.89 , 10 and 100 kg(m -3 s -1 ) revealed that the center of mass of the thermal dose due to focal volume heating of each simulated sonication appeared to be in approximately the same location regardless of the value of perfusion, while the values of the focal volume thermal dose distribution became closer in magnitude to that of the nearfield thermal dose distribution as perfusion was increased (Supplementary Figure S2). Regarding these five sonications, when comparing the location of the center of mass of a thermal dose distribution of the sagittal slice of a sonication simulated using a value of perfusion of 1.89 kg(m -3 s -1 ) to that of a higher value the resulting displacement values were large at least for the highest value of omega (see previous paragraph). These observations support the aforementioned hypothesis that for a sonication, a case where simulated focal volume heating is low enough to yield a focal volume thermal dose distribution comparable in magnitude to that of the nearfield, can yield a sagittal slice thermal dose distribution center of mass showing a large difference when compared to the center of mass of a simulated sagittal slice thermal dose distribution showing high focal volume heating relative to the nearfield while the locations of the region containing locally high values of thermal dose due to focal heating are approximately the same in both simulated sagittal slice distributions. This effect could be further analysed by recalculating the center of mass in a cropped section of the sagittal simulated MRthermometry slice instead of the full simulated slice, in which case the aforementioned nearfield heating contribution to the shift would likely be greatly reduced.   Figure S 4 The idea that intra-patient variation of tissue parameters may be a practically significant contribution to the ranges of discrepancies between peak temperature curves of simulations and treatment is further supported by performing a comparison of the power normalised initial peak temperature heating rates of simulations and treatment. For example, within some patients there are sonications having treatment cell positions with the same depth coordinate that exhibit a range of patient treatment peak initial heating rate results that is roughly equal to the range of analogous results of simulations when varying alpha between 1.2 and 8.6 Np(m -1 MHz -1 ). The simulated heating rates in this figure were estimated using the difference of the peak temperature of the second simulated sagittal MR-thermometry acquisition and the baseline temperature of 37 degrees Celsius. This figure also illustrates the low degree of clinical heating exhibited by patient 3. In addition, of the 67 sonications considered in this simulation study, the power normalised initial heating rate estimated using the peak temperature in the focal volume of the first and third sagittal MR-thermometry slice acquisitions averaged over sonications was on average lower for patient 3 compared to the other patients, despite the sonications of patient 3 having treatment cell position depth coordinates that were similar to as well as anterior of those of sonications exhibiting a higher degree of heating based on the aforementioned metric that were performed in the treatment of another patient.  2,4.9, MHz -1 ) respectively where the circle, square, and triangle denote perfusion values of 1.89, 10, and 100 kg(m -3 s -1 ) respectively. The indexing of sonications in this figure is independent of the indexing elsewhere in the manuscript and supplementary section. The perfusion parameter value of 100 kg(m -3 s -1 ) might be too high to be representative of the patient perfusion in many of the simulated sonications in which it was utilised since the peak temperature curves of the simulations utilizing a fibroid absorption value of 8.6 Np(m -1 MHz -1 ) and this value as a perfusion parameter typically exhibited a high degree of cooling larger than that of the peak temperature curves of patient treatment. Sonication 8 of patient 2 demonstrates good agreement with a perfusion parameter of 10 kg(m -3 s -1 ) where varying the perfusion parameter in simulations from 1.89 to 10 kg(m -3 s -1 ) did not have a strong effect on peak temperature curves, but the effects on the 240EM43 thermal dose threshold volumes are substantial.  There is data for 19 sonications in the bottom panel because one of the six sonications of patient 2 had a simulated 240EM43 thermal dose threshold volume of zero for both of the aforementioned values of perfusion, thus yielding an undefined ratio. This figure shows that while varying the perfusion parameter from 1,89 to 10 kg(m -3 s -1 ) does not have a large effect on peak temperature curves, the effect on 240EM43 thermal dose threshold volumes can still be substantial. Figure S 9: The coordinate wise differences between the simulated and measured thermal dose center of mass are shown for N=67 sonications using the solid blue narrow boxplots with outliers shown using blue circles. For example, a positive value for an AP displacement means that the simulation of a sonication resulted in a thermal dose center of mass posterior of the treatment. The corresponding results of the analysis in which the simulated sagittal MR-thermometry slice and simulated stack of coronal MR-thermometry slices was placed at the treatment cell position with the misregistration correction added to it are shown using the wide hollow box plots with outliers shown using red plus signs. The value of absorption used is denoted by α in units of Np(m -1 MHz -1 ) and the value of perfusion used was ω = 1.89 kg(m -3 s -1 ). Figure S 10: Ratio of sagittal MR-thermometry peak temperature rise of simulation to patient treatment as a function of time averaged over sonications within individual patient treatments as well as throughout patient treatments is shown using solid lines. The blue, black, and red colors each represent a value of fibroid absorption coefficient of 1.2, 4.9, and 8.6 Np(m -1 MHz -1 ) respectively that was used in simulations. The data is shown only for the durations of the innermost trajectory in order for the averaging to be meaningful (for example averaging a sonication that is cooling while another is heating at a particular timepoint may not be particularly meaningful). The number of sonications is thus variable as a function of time because the sonications do not all have the same innermost trajectory duration. To indicate the degree of discrepancy between simulation and treatment as well as information about the distribution of the data, notched boxplots in solid colors are shown for each timepoint. The triangles represent the notches of the boxplot calculated using the formula q2±1.57 (q3-q1)/(N 0.5 ) where q1, q2, and q3 are the 25th, 50th, and 75th percentiles respectively and N is the number of observations. The non-overlapping intervals indicated using triangles indicate statistically significant differences of the medians at 5% significance level. The corresponding results of the analysis in which the simulated sagittal MR-thermometry slice was placed at the treatment cell position with the misregistration correction added to it are shown using dotted lines and hollow boxplots. The horizontal constant dashed line having black color denotes unity. Np(m -1 MHz -1 ). Bottom: A Bland-Altman plot with ratios of 240EM43 volumes of simulations to treatment as a function of the average of the two quantities of N=54 sonications each having a patient treatment 240EM43 volume greater than 2 voxels. The limits of agreement are mean ratio ± 2 standard deviations.